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Abstract  Our  effort  has  been  continued  for  developing  and  analyzing  high  order  methods 
for  the  interface  problems.  Specifically,  We  have  developed  an  adaptive  mesh  refinement 
(AMR)  for  the  Immersed  Boundary  (IB)  and  the  Immersed  Interface  Method  (IIM)  using 
the  level  set  method  for  elliptic  interface  problems,  the  sensitivity  calculation  with  respect 
to  the  parameter  and  shape  for  flow  past  obstacle  and  multi-moment  immersed  interface 
method  for  wave  equation  and  Hamilton- Jacobi  equation.  Our  approach  uses  the  physical 
augmented  variables  and  the  domain  embedding  technique  and  allows  us  to  develop  fast  and 
robust  implementation  of  the  immersed  interface  method  for  general  interface  problems.  In 
particular,  the  sensitivity  calculation  with  respect  to  the  Reynolds  number  for  the  vortex 
shading  in  the  flow  past  obstacle  problem  near  the  critical  Reynolds  number  (Re=47)  shows 
the  enhanced  and  detailed  sensitivity  toward  the  instability.  The  sensitivity  calculation 
provides  more  precise  and  detailed  sensitivity  of  the  solution  and  describes  the  dynamical 
change  due  to  the  variation  in  the  Reynolds  number.  The  immersed  interface  method 
accurately  computes  the  velocities  and  the  stress  on  the  interface. 

This  is  the  first  time  that  an  AMR  has  been  applied  to  the  immersed  interface  method. 
The  new  algorithm  is  based  on  the  level  set  representation  of  the  interface.  This  unique  fea¬ 
ture  enables  us  to  refine  the  mesh  only  in  the  neighborhood  of  the  interface.  The  results  for 
elliptic  interface  problems  with  arbitrary  shapes  are  tested.  An  accurate  sensitivity  analysis 
with  respect  to  the  change  in  the  interface  condition  and  geometry  has  been  computed  effi¬ 
ciently  via  the  immersed  interface  method  as  the  forward  and  the  sensitivity  equation  solver. 
We  develop  and  analyze  a  new  multi-moment  method  for  one-dimensional  hyperbolic  equa¬ 
tions  with  discontinuous  coefficients.  The  method  is  based  on  the  backward  characteristic 
method  and  uses  the  solution  and  its  derivative  as  unknowns  and  cubic  Hermite  interpola¬ 
tion  for  each  computational  cell.  An  exact  update  formula  for  solution  and  its  derivative  for 
variable  wave  speed  is  derived  and  used  for  efficient  time  integration.  At  points  of  discon¬ 
tinuity  we  develop  a  piecewise  cubic  Hermite  interpolation  based  on  interface  conditions. 
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The  method  is  extended  to  the  one-dimensional  Maxwell  equation  with  variable  material 
properties.  The  extension  of  methods  to  multi-dimensional  case  is  under  progress,  in  which 
we  use  Poisson  formula  for  wave  equation  as  a  building  block. 

1  Project  Summary 

Our  effort  has  been  continued  for  developing  and  analyzing  high  order  methods  for  the  in¬ 
terface  problems.  We  have  developed  an  adaptive  mesh  refinement  (AMR)  for  the  Immersed 
Boundary  (IB)  and  the  Immersed  Interface  Method  (IIM)  using  the  level  set  method  for 
elliptic  interface  problems.  We  have  continued  the  sensitivity  calculation  with  respect  to  the 
Reynolds  number  for  the  vortex  shading  in  the  flow  past  obstacle  problems  and  a  accurate 
sensitivity  analysis  with  respect  to  the  change  in  the  interface  condition  and  geometry  has 
been  computed  efficiently  via  the  immersed  interface  method  as  the  forward  and  the  sensi¬ 
tivity  equation  solver.  We  have  developed  and  analyzed  a  new  multi-moment  method  for 
one-dimensional  hyperbolic  equations  with  discontinuous  coefficients.  The  method  is  based 
on  the  backward  characteristic  method  and  uses  the  solution  and  its  derivative  as  unknowns 
and  cubic  Hermite  interpolation  for  each  computational  cell.  An  exact  update  formula  for 
solution  and  its  derivative  for  variable  wave  speed  is  derived  and  used  for  efficient  time  inte¬ 
gration.  At  points  of  discontinuity  we  develop  a  piecewise  cubic  Hermite  interpolation  based 
on  interface  conditions.  The  method  has  been  applied  to  the  Hamolton-Jacobi  equation. 

2  Status/Progress 

Elliptic  interface  problems  are  fundamental  problems  for  incompressible  flow  solvers  for 
Stokes  flow  and  Navier-Stokes  equations.  This  is  the  first  time  that  an  AMR  has  been 
applied  to  the  immersed  interface  method.  The  new  algorithm  is  based  on  the  level  set 
representation  of  the  interface.  This  unique  feature  enables  us  to  refine  the  mesh  only  in 
the  neighborhood  of  the  interface.  The  results  for  elliptic  interface  problems  with  arbitrary 
shapes  are  excellent.  Let  ip(x)  be  a  levelset  function  whose  zero  level  set  {ip  =  0)  defines 
the  interface.  If  we  have  a  uniform  Cartesian  grid  with  mesh  size  h\.  For  the  grid  points 
(: Xi,yj )  in  the  tube  \(p\  <  5 ,  we  generate  a  finer  grid  with  smaller  mesh  size  ( h-2  <  hi). 
Generate  the  finite  difference  equations  using  the  standard  5-point  or  9-point  stencil.  The 
procedure  can  be  repeated  hierarchically.  Whenever  the  finite  difference  scheme  at  a  grid 
point  in  the  new  mesh  needs  values  at  points  that  are  not  grid  points  in  any  level  of  the 
mesh  (hanging  nodes),  an  interpolation  scheme  is  used.  The  resulting  linear  system  of  finite 
difference  equations  is  solved  by  an  algebraic  multi-grid  solver. 

The  AMR  features  include: 

•  Better  accuracy  near  the  interface 
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Table  2.1:  Comparison  of  results  using  a  uniform  mesh  and  that  of  the  AMR-IIM  with 
the  same  finest  resolution  h.  The  runtimes  in  second  by  algebraic  multigrid  solver  are  also 
provided. 


Uniform 

AMR(+1)[4] 

m 

# 

Error 

Ratio 

CPU 

m 

A 

# 

Error 

Ratio 

CPU 

40 

1521 

8.44e-4 

.03 

20 

1 

805 

9.07e-4 

.02 

80 

6241 

2.45e-4 

3.4 

.08 

40 

2 

2745 

2.86e-4 

3.2 

.06 

160 

25281 

6.69e-5 

3.7 

.31 

80 

4 

10201 

7.80e-5 

3.7 

.19 

320 

101761 

1.57e-5 

4.3 

1.3 

160 

8 

39225 

1.89e-5 

4.1 

.67 

640 

408321 

3.65e-6 

4.3 

5.4 

320 

16 

153805 

5.04e-6 

3.8 

2.5 

•  Less  degree  of  the  freedom,  i.e.  smaller  system  to  solve 

•  Combined  with  the  level  set  method  so  it  is  flexible  in  dealing  with  complicated  domain 

•  Use  an  algebraic  multi-grid  solver  which  is  fast. 

We  have  applied  to  an  elliptic  interface  problem: 

uxx  +  uyy  =  0,  on  0\T,  [u]  r  =  0, 

Figure  2.1  shows  the  three  level  refinement  for  a  eight  star  shape  and  Table  2.1  shows  the 
AMR  performance  for  the  corresponding  solutions. 
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Figure  2.1:  A  two  level  AMR  for  a  eight  star  shape 

2.1  The  sensitivity  analysis  and  flow  control 

We  haveinvestgated  the  sensitivity  analysis  with  respect  to  the  Reynolds  number  for  the 
flow  past  obstacle  problem.  To  carry  out  such  analysis,  at  each  time  step,  we  need  to  solve 
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the  incompressible  Navier-Stokes  equations  on  irregular  domains  twice,  one  for  the  primary 
variables;  the  other  is  for  the  sensitivity  variables  with  homogeneous  boundary  conditions. 
The  Navier-Stokes  solver  is  the  augmented  immersed  interface  method  for  Navier-Stokes 
equations  on  irregular  domains  which  is  developed  under  the  current  grant.  Our  sensitivity 
analysis  can  predict  the  critical  Reynolds  number  at  which  the  vortex  shading  begins  to 
develop  in  the  wake  of  the  obstacle.  Numerical  experiments  are  shown  to  illustrate  how  the 
critical  Reynolds  number  varies  with  different  geometric  settings. 

The  sensitivity  equation  with  respect  to  the  constant  viscosity  y  for  the  incompressible 
Navier  Stokes  equation  is  given  by;  v  =  ^u,  q  =  ^ p  satisfy 


,<9v  _  _  .  _  ^  g(Vv  +  VvT) 

p(—  +  u  •  Vv  +  v  •  Vu)  +  \7q  =  V  - - - 


+  V- 


(Vu  +  Vu5 


V  •  V  =  0. 

We  have  carried  out  our  sensitivity  analysis  to  examine  the  effect  of  the  shape  and 
orientation  of  the  obstacle  on  the  vortex  shading.  The  obstacle  that  we  considered  is  an 
ellipse  whose  level  set  function  is 

<p(x,  y )  =  s/(x  +  7)2/0.62  +  y2/0-42  -  1 

with  different  orientations.  We  see  that  critical  Reynolds  number  of  Re*  =  l/y*  is  quite 
different  for  different  orientations.  The  sensitivity  analysis  does  give  a  good  estimate  of  the 
critical  number  at  which  the  vortex  shading  is  about  to  develop.  In  Figure  2.2  (left),  we 
show  the  sensitivity  analysis  of  an  ellipse  at  different  orientations.  The  plots  on  the  left 
are  the  plots  of  the  vorticity  (top)  and  sensitivity  vorticity  (bottom)  strengths  versus  the 
parameter  1/y  ~  Re.  From  the  plot,  we  can  predict  that  the  critical  Reynolds  number  is 
significant  different  with  the  orientation.  The  results  are  also  agree  with  intuition.  With 
the  orientation  is  parallel  to  the  flow  direction,  the  flow  in  the  wake  of  the  obstacle  is  more 
stable  than  that  against  the  flow  direction.  As  we  can  expected,  we  can  predict  that  the 
critical  Reynolds  number  is  the  largest  for  the  horizontal  case  (around  70),  decreases  for  the 
oblique  ellipse  (around  45),  and  the  smallest  for  the  vertical  case  (around  40).  The  plots 
on  the  right  are  vorticity  (top)  and  sensitivity  vorticity  (bottom)  at  some  numbers  which 
supports  our  discussions  above. 

We  have  also  carried  out  the  sensitivity  analysis  for  a  flow  past  two  separated  cylinders. 
The  level  set  function  is  defined  as  follows 

Vi(x,  y)  =  y/(x  +  7)2  +  y2  -  1,  <p2(x,  y)  =  y/{x  -  x2)2  +  y2  -  1, 

<pi(x,y)  =mm{tp1(x,y),(p2(x,y)  }. 

In  Figure  2.2  (right),  we  show  the  sensitivity  analysis  for  the  two  stationary  cylinders.  We 
show  three  different  cases.  In  the  first  case  (a),  the  two  cylinders  are  centered  at  (—7,  0) 
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(a),  Horizontal  ellipse,  the  right  plot  are  vorticity  or  saisitivity  vortidty  at  Re  =  70. 


(b),  Oblique  ellipse,  the  right  plot  are  vorticity  or  sensitivity  vorticity  at  Re  =  45. 


€Bgsy 

(c),  Vatical  ellipse,  the  right  plot  are  vorticity  or  sensitivity  vorticity  at  B/t  =  40. 


(a),  Two  cylindas  centered  at  (-7,  0)  and  (-4,  0).  The  critical  Reynolds  is  about  54. 


(b),  Two  cylinders  catered  at  (—7,  0)  and  (—3,  0).  The  critical  Reynolds  is  about  64. 


(c),  Two  cylinders  centred  at  (-7,  0)  and  (0,  0).  The  critical  Reynolds  is  about  65. 


Figure  2.2:  Vorticity  strength  versus  the  Reynolds:  The  top  row  is  for  the  primitive  and 
the  bottom  row  is  for  the  sensitivity  equations. 


and  (—4,  0),  respectively.  They  are  considered  as  relatively  close.  We  see  that  the  critical 
Reynolds  number  is  about  Re*  =  54  which  is  larger  than  that  of  a  single  cylinder.  As  we 
increase  the  distance  between  two  cylinders  by  one  unit,  we  see  that  the  critical  Reynolds 
number  increases  by  10,  see  Fig.  ??  (b)  for  which  the  two  cylinders  are  centered  at  (—7,  0) 
and  (—3,  0)  respectively.  If  we  further  increase  the  distance  between  two  cylinders  by 
another  3  units,  the  critical  Reynolds  number  remains  almost  the  same.  In  Fig.  ??,  the  left 
plots  are  the  strength  of  the  vorticity  and  sensitivity  versus  the  Reynolds  number  defined 
as  1  hi.  The  right  plots  are  the  vorticity  and  sensitivity  plot  at  t  =  100  and  a  Reynolds 
number  that  is  close  to  the  critical  Reynolds  number  54,  64,  and  65  respectively. 

We  have  investigated  with  rotating  objects  for  the  single  obstacle  and  two  separated 
cylinders.  Our  sensitivity  analysis  can  provide  a  good  estimate  of  the  critical  number  at 
which  the  vortex  shading  is  about  to  develop  even  for  dynamical  flow  cases.  We  will  continue 
to  investigated  for  the  different  flow  regime  (the  both  backward  and  forward  steps  and  cavity 
flows)  under  various  applied  flow  controls  and  develop  the  sensitivity  analysis  tool  for  the 
actuator  and  sensor  dynamics. 
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2.2  Multi-moment  method  for  hyperbolic  equations 

We  have  developed  a  multi-moment  method  for  the  wave  propagation  in  discontinuous 
media.  We  considered  and  transport  as  well  as  advection  equation  as  a  test  model: 

ut  +  c{x)ux  =  0,  ut  +  ( c(x)u)x  =  0,  x  £  R.  (2-1) 

where  the  wave  speed  c  =  c{x)  >  0  is  variable  and  discontinuous.  Our  method  is  motivated 
and  closely  related  to  CIP  methods  [?] .  CIP  is  one  of  the  numerical  methods  that  provides  an 
accurate,  less-dispersive  and  less-dissipative  numerical  solution.  Our  method  uses  the  exact 
integration  in  time  by  the  characteristic  method  and  uses  the  cubic  Hermite  interpolation 
in  each  cell  [xj-i,Xj\  based  on  solution  values  and  its  derivatives  at  two  endpoints  Xj-i,  Xj. 
That  is,  we  develop  the  exact  solution  formula  for  solution  and  its  derivative  for  (2.1).  The 
Hermite  cubic  interpolation  is  then  used  to  evaluate  the  formula  locally.  Our  development 
of  CIP  is  entirely  based  on  the  characteristic  method  and  results  in  a  different  (improved 
and  simpler)  scheme  than  the  conventional  one  for  equations  with  variable  wave  speed. 
For  example  the  new  method  allows  us  to  take  an  arbitrary  time  step  (no  CFL  limitation) 
without  losing  the  stability  and  accuracy.  We  analyzed  the  von-Neumann  stability  of  the 
proposed  method  for  a  constant  speed.  We  develop  an  immersed  interface  method  [?]  for  the 
discontinuous  wave  speed  case,  i.e. ,  we  construct  a  piecewise  cubic  Hermite  interpolation  at  a 
cell  [xj-i,Xj\  which  contains  a  point  of  discontinuity  of  c(x)  using  proper  interface  conditions 
and  solution  values  and  its  derivatives  at  two  endpoints  Xj- 1,  Xj.  This  interface  treatment 
is  applied  for  the  one-dimensional  Maxwell  equation  with  variable  material  properties.  We 
first  approximate  the  variable  coefficients  by  a  piecewise  constant  (discontinuous)  coefficient. 
The  d’Alembert’s  based  method  for  the  Maxwell  equation  that  extends  our  characteristic 
based  method  to  Maxwell  system  is  developed  for  the  piecewise  constant  media  and  then 
applied  to  Maxwell  system  with  piecewise  constant  coefficients. 

In  the  case  of  the  transport  equation  we  have 

c(y ) 

u(t  +  At,  x)  =  u(t,y),  ux  (t  +  A  t,x)  =  ——ux(t,,y) 

c(x) 

where  y  is  the  backward  characteristic  of  x  €  R  determined  by 

— x(t)  =  c(x(t)),  x{t  +  8t)  =  x  and  y  =  x(t). 

Suppose  c(x)  is  piecewise  constant  and  discontinuous  at  x*  £  (xj-\,Xj),  We  construct  the 
piecewise  cubic  interpolation  F±(x )  at  each  time  step  tn: 

F±(x )  =  ^2ak(x~x*)k 

k= 0 
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Figure  2.3:  Discontinuous  wave  speed:  advection  case 

in  each  interval  [xj- i,x*]  and  \x*,Xj\,  respectively.  We  determine  the  eight  unknowns  via 
the  interface  relations  [it]  =  [cux]  =  [ c2uxx ]  =  [c3uxxx\  =  0  and  and  the  interpolation 
conditions  u(xj-i)  =  ix"_ l5  , ux(xj-\ )  =  i>”_,  at  Xj-\  and:  u(xj)  =  it",  ,ux(xj)  =  v™  at  Xj. 
Thus,  the  update  formula  becomes 

un+l  =  p+  ^  _C+Af^  yn+ 1  =  £_  F+  (Xn  _  c+  At ) . 

In  Figure  2.3  solutions  in  discontinuous  wave  speed  media  for  the  advection  equation. 

We  are  currently  extending  our  treatment  to  the  multi-dimensional  wave  equations.  The 
building  block  of  the  extension  is  the  Poisson  formula  for  solutions  to  the  wave  equation, 
i.e.  one  can  derive  an  exact  time  integration  method  in  time  for  the  homogeneous  media 
case.  Also,  we  uses  the  (cell-wise)  bi-cubic  interpolation  based  on  moments  u,  ux,  uy,  uxy. 
We  analyze  the  von-Neumann  stability  of  the  proposed  method  and  we  establish  CFL 
number  is  one  for  the  proposed  method.  Despite  its  complexity,  the  method  offers  a  highly 
efficient  method  to  attain  the  desired  accuracy  due  to  the  relaxed  CFL  number  limitation 
and  accurate  space  resolution  by  the  bi-cubic  profile.  The  method  preserves  discontinuous 
profiles  very  accurately  without  any  smearing  and  distortion,  e.g.  see  Figure  2.4.  Also,  the 
method  computes  directly  the  physical  quantities,  e.g.,  current  and  electric  field  gradient, 
very  accurately.  The  immersed  interface  treatment  will  be  integrated  to  develop  a  higher 
order  fully  time-space  method  for  heterogeneous  media. 

2.3  Multi-moment  method  for  Hamilton-Jacobi  equation 

We  develop  a  multi-moment  approximation  method  for  Hamilton  Jacobi  (HJ)  equations, 
which  arise  in  many  applications  such  as  geometrical  optics,  crystal  growth,  etching,  com- 
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Figure  2.4:  Scattering  from  a  square  (CFL=1) 


puter  vision,  obstacle  navigation,  path  planning,  photolithography,  and  seismology  and  con¬ 
trol  problems.  In  general,  the  solutions  usually  develop  singularities  in  their  gradient  even 
with  smooth  initial  conditions.  In  many  application  of  HJ  it  is  very  essential  to  compute  the 
gradient  of  the  solution  along  with  solution  accurately.  Especially,  for  control  problems  the 
optimal  feedback  solution  is  determined  using  the  gradient  of  the  value  function  (solution) . 

We  develop  approximation  methods  based  on  the  Lax-Hoph  formula  for  the  convex 
Hamiltonian  and  the  maximum  principle  for  the  control  problem  as  the  (exact)  time  in¬ 
tegration.  The  formula  calculates  the  unique  viscosity  solution  of  HJ  equations.  In  order 
to  update  the  solution  to  the  next  time  step  we  use  a  locally  defined  quadratic  and  cubic 
interpolation  of  the  solution  based  on  the  solution  value  and  its  gradient  at  each  square 
local  cell.  To  complete  the  multi-moment  method  we  develop  the  update  formula  for  the 
gradient. 

For  Hamilton-Jacobi  equation 


Vt  +  H(vx)  =  0,  i£H, 

where  FI  is  a  convex  Hamiltonian,  for  example  H(P )  =  —  \p\ 2 ,  \p\.  Let  H*  be  the  convex 
conjugate  of  H\ 

H*(v)  =  sup  ((v,p)  -  H(p)). 

p 

We  have  the  solution  update  formula: 

V(z,t  +  At)  =  minCen  {V((,t)  +  AtH*(^f)}. 


VV{z,t  +  At)  €  dH*  (*=&■). 
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Figure  2.5:  ID  Hamilton-Jacobi  equation:  u  =  VX 


where  £ x  is  a  minimizer. 

The  following  is  our  proposed  algorithm: 

•  Initialize 

•  Solve  the  two-sided  quadratic  minimizations 

Ui  =  argmin {Fi-lti(tn,y)  +  A tH*(^r)}  over  ;*] 

and 

yf  =  &rgmm{Fiji+1(tn,y)  +  A tH*{B^L)}  over  [xi,xi+1] 

•  Find  y f  (either  y^~  or  y^~)  that  achieves  the  minimum  of  the  two-sided  minimization; 

vi+1  =  mi  n(iJ,t-i,i(in,yr).-FM+i(*n>yi')) 


l  t  _  (ll'l 

..ri+1  _  U _ 

“  At 


Here  Fl_\l  and  Ftj+\  are  the  interpolation  function  on  each  interval  [xi-\,Xj\  and 
\xi,Xi+ 1],  respectively.  The  building  block  of  our  integration  method  for  Hamilton-Jacobi 
equation  is  the  Hoph-Lax  formula  and  variant  of  methods  can  be  developed  based  on  var¬ 
ious  proper  interpolation  method  for  the  local  solution  profile.  We  have  also  successfully 
applied  the  method  for  Hamilton-Jacobi  equation  associated  with  optimal  control  problems. 
In  Figure  2.5  solutions  are  shown  for  the  case  H(p)  =  for  square  wave  u  =  Vx  (Burgers 
equation) . 
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